On the trapping reaction with mobile traps 
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We present Monte Carlo results for the two-species trapping reaction A + B — > B with diffusing 
A and B on lattices in one, two and three dimension. We use a novel algorithm which permits to 
simulate survival probabilities of A particles down to < 10 -30 with high accuracy. The results for 
the survival probability agree much better with the exact asymptotic predictions of Bramson and 
Lebowitz (Phys. Rev. Lett. 61, 2397 (1988)) than with the heuristics of Kang and Redner (J. Phys. 
A 17, L451 (1984)). But there are very large deviations from either which show that even these 
simulations are far from asymptotia. This is supported by the rms. displacement of A particles 
which clearly show that the asymptotic regime has not been reached, at least for d = 2 and d — 3. 
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Even the simplest reaction diffusion systems with only 
two-body reactions and without particle production show 
very rich and not yet fully understood phenomena. One 
prototype is the absorption of diffusing particles on ran- 
domly located static sinks. This model, studied first by 
Smoluchowski nearly hundred years ago , shows an ex- 
ponential decrease of the particle density in mean field 
theory, but the exact asymptotic behaviour found by 
Donsker and Varadhan M is a stretched exponential in 
any finite dimension d. Another classic example, the re- 
combination A + A — > A, was studied by Feller ||. It 
leads to tia ~ <~ 1//2 in d = 1, as compared to the mean 
field solution ua ~ t^ 1 which holds for d > 2. 

Particle annihilation with two mobile species A and 
B, either of the type A + B^B or A + B^O, were 
studied by Ovchinikov and Zeldovich E] and Toussaint 
and Wilczek ||. The motivation of stemmed from 
the fate of magnetic monopoles in the early universe, but 
applications to condensed matter physics and chemical 
kinetics are of course more numerous. 

If one starts out with equal densities for A and B, 
n ^i(0) = ns(0), then the reaction A + B — > leads to 
n ~ i -1 / 4 ||. This is different from the cases ua(0) < 
nfl(O) and A + B — > B. In both these cases ua <C ub for 
late times, and one expects identical asymptotics. This 
was indeed proven rigorously by Bramson and Lebowitz 
Ipl who obtained 



n A (t) 



exp(— Ai \/t) 
exp(— Xzt/ lni) 
exp(-A d i) 



d = 1 
d= 2 
d > 3 



(1) 



with unknown constants A^. In contrast to this, different 
asymptotics had been predicted for A + B — > B and for 
A + B — > with Ua(0) < ng(Q) by Kang and Redner 
jj], ^ by heuristic arguments. 

Verifying Eq.(|l|) numerically has turned out to be 
about as difficult as verifying the Donsker- Varadhan 
stretched exponential. The first simulations by Kang and 
Redner |?j agreed with their own (supposedly wrong) 
heuristic asymptotics. Subsequent simulations [Hj 
were judged even by the authors as inconclusive. The 
main problem seems to be that there are large finite 
time corrections to the asymptotic behaviour. Thus one 



would like to simulate up to very long times. But with 
the straightforward approaches used so far it is practi- 
cally impossible to estimate survival probabilities smaller 
than ss 10 -8 , even with the most powerful present day 
computers. 

It is the purpose of this work to present simulations for 
the reaction A + B — * B which go far beyond this. We 
shall not give any results for A + B — * since our special 
numerical methods cannot be applied to that case. More 
specifically we consider regular lattices on which the B 
particles perform independent random walks. A particles 
also perform random walks until they hit upon a B par- 
ticle in which case they are instantaneously absorbed. 
Initial conditions are such that all Bs are uncorrelated 
with homogeneous concentration c. Since A particles do 
not interact with each other, their concentration is ir- 
relevant. In the actual simulations we shall use cither 
one or two As in the initial configuration. Notice that 
in this model the B particles act as catalyzers, and their 
distribution is Poissonian at any time, if they start out 
independently at t — 0. 

Our algorithm is related to an algorithm used re- 
cently |Tl| for the trapping (Donsker- Varadhan) prob- 
lem. There we were able to clarify the cross-over from the 
mean field type to the stretched exponential behaviour. 
In the present case we shall see that - despite going to 
much longer times than in earlier simulations - we still do 
not yet fully understand the cross-over to the Bramson- 
Lebowitz asymptotics. 

The algorithm has several essential ingredients. The 
first is that we use cloning ( "enrichment" ) of configura- 
tions with surviving A particles and a bias for the diffu- 
sion of A such that less of them are absorbed. This bias 
is compensated by suitably chosen weights, i.e. we al- 
ways deal with nontrivially weighted ensembles. On the 
other hand, the simulation is stopped as soon as one of 
the As is absorbed, or if the weight of the configuration is 
too small. Thus all our results are based on conditional 
probability distributions, conditioned on the survival of 
all As. These features are implemented by PERM 
which is a general growth method (using 'sequential im- 
portance sampling with resampling' in the sense of | fl3"| ) 
and has been very successful in a large number of prob- 
lems ini pi. 
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Assume that we have a single A (the case for k A par- 
ticles with k > 1 is straightforward) which has arrived 
at site i at time t. The fact that we condition on those 
events where A is not absorbed means that there cannot 
be a B at this site. Thus the homogeneity of the uncon- 
ditioned distribution is broken, and an effective hole in 
the B distribution is introduced. For times > t this hole 
makes a random walk. If it meets a hole produced at a 
time 7^ t, the two recombine. 

In principle one could simulate these holes explicitly, 
i.e. one could simulate the As in a background of B 
where a B particle (sic!) is removed each time it hits an 
A. We indeed did perform such simulations. They agreed 
with straightforward simulations not using PERM or any 
conditioning and were more accurate, but the accuracy of 
both type of simulations was rather poor in comparison 
with our final algorithm. 

Our final ingredient is that we do an exact summation 
over the B paths. This depends crucially on the fact that 
the B distribution is Poissonian. The latter is still true 
for the conditioned distribution (it would not be true in 
the reaction A + B — ► 0, therefore our method cannot be 
used for it). A Poisson distribution is uniquely charac- 
terized by its mean. The evolution of this mean is de- 
scribed by a modified diffusion equation. More precisely, 
we describe the B density by a homogeneous background 
c minus a density p(i,t) of holes. The latter is at every 
time step set equal to c at the actual A position, but oth- 
erwise evolves according to the simple diffusion equation 
p(i, t + 1) = (2d) -1 Yl<j i> p{ii I ts initial condition is 
p(i, 0) = 0. At each time step the A (which is assumed to 
be at site i) has a chance exp(p(i, i) — c) to be absorbed, 
i.e. the weight of the event in the conditioned ensemble 
decreases by a factor exp(c — p(i, t)). The algorithm used 
in |H| for the trapping problem is just a simplification 
where we omit the diffusion of the holes. 

All our simulations were done on workstations, with 
a total of a few hundred hours of CPU time. All results 
were carefully checked for small t against straightforward 
brute force simulations with the most simple algorithm. 
In addition we also made simulations with algorithms of 
intermediate complexity and efficiency 

We first discuss results for d = 1. Here we can use 
lattices so large that none of the holes ever reaches the 
boundary, thus we have no finite size effects at all. In 
Fig. la we show the survival chances Pa{^) of a single A 
particle and of a pair of particles which started at the 
same site. Although they become as small as 10 -35 , the 
statistical errors are much smaller than the thickness of 
the lines. The fact that Pa for a pair is larger that the 
square of Pa for a single particle is easily understood: 
If already one A particle has survived in some region, 
there are less than average B particles in this region, 
and the second A particle has a bigger survival chance. 
Theoretically || we expect Pa(£) ~ exp(— const y/i), 
at least for a single particle. In Fig. lb we thus show 
VtlnPyi(t) on a semi- logarithmic scale. Error bars are 
here < 0.001. Thus the decline for t > 1000 is statis- 
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FIG. 1: (a) Semi-logarithmic plot of the survival probability 
in d = 1 for c = 0.5. The upper curve is for a single A particle, 
the lower for a pair. 

(b) — \ft In Pa (t) from the same data as in panel (a), plotted 
on a logarithmic horizontal scale. 



tically highly significant for both curves. Unless we ac- 
cept that there is an error in H, we have to conclude 
that at least this decline of — vfln P^(t) for a single A 
particle is a finite t effect, and that the true asymptotic 
behaviour sets in at t ^> 10 4 . We also compared our data 
with the prediction — lnP^(f) ~ i -1 / 4 of 0, with even 
worse agreement. Although that latter prediction was for 
A + B — > B, ua{0) < ns(0), we pointed out already that 
both models should show the same asymptotics. We thus 
conclude that the data are in rough agreement with 
but the very big deviations are surprising (in particular 
since they are not monotonic!) and not understood. 

The mean squared displacement (R\) of the A par- 
ticle is shown in Fig. 2. In this figure we also show the 
squared distance between two A particles which started 
off simultaneously at the same site. We see that AR in- 
creases less fast than R, indicating an effective attraction 
mediated by the B particles: If already one A particle has 
survivied in some region, there are less than average B 
particles in this region, and the second A particle will 
not only survive longer, but it will also survive preferen- 
tially in a region close to the first one. But this appears 
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FIG. 2: Full line: mean squared displacement of the A particle FIG. 4: Log-log plot of mean squared displacements in 1, 2 
in d — 1, for c = 0.5. Broken line: mean squared distance be- and 3 dimensions, divided by y/t. 
tween two A particles, when conditioned on survival of both. 
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FIG. 3: Full line: -t _1 \nP A {t) for d = 2 and c = 0.5; Broken 
line: — i -1 hit lnPA(t) from the same data. The horizontal 
scale is logarithmic, the vertical is linear. 



to be a weak effect. The most economic conclusion from 
Fig. 2 is that both curves converge to the same scaling be- 
haviour, R 2 ~ (Ai?) 2 ~ t v with v = 0.5 to 0.6. But any 
determination of a critical index should be taken with 
very big caution in view of Fig.l: we should not take any 
behaviour seen for t < 10 5 as reflecting the asymptotic 
behaviour. 

With increasing d, updating p(i,t) becomes more and 
more time consuming. Thus we can simulate only for 
much shorter times, if we want to avoid excessive CPU 
times or large finite size corrections. For d = 2, results 
for i -1 lnPA(t) with c = 0.5 are shown in Fig. 3, together 
with — t^ 1 hit lnP^i) . For large values of t the data 
agree much better with Eq.|l]) than with the alternative 
prediction - \nP A (t) ~ at-bt 1 / 2 of @. The factor 1/lnt 
in Eq.(|l|) seems to be correct asymptotically, although it 
makes agreement worse for small t. Anyhow, deviations 
from Eq.([|) are substantial, and even for the largest t 
where the curve appears to be horizontal in the figure 
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FIG. 5: Log-linear plot of -r 1 P A (t) for d = 3, c = 0.5. 



(and where, by the way, P A {t) ~ 10~ 98 ) it still shows a 
definite curvature. 

The fact that asymptotia cannot have been reached by 
these 2-d data is most clearly seen from (R A ) ■ It is shown 
in Fig. 4, together with the analogous data from d = 1 
and d = 3. To make the point particularly clear we show 
there (R\)/\/i against t. If d = 2 is the upper critical 
dimension for A + B — * B as suggested by Eq. ([!]), then 
we should expect (R\) ~ t up to logarithmic corrections 
in d = 2. We also should expect that the data for d = 2 
should fall between those for d = 3 (where we expect 
(R 2 A ) ~ t) and d = 1 (where (R 2 A ) - t v with v < 1). 
Figure 4 shows a completely different behaviour. On the 
one hand, (R\) is very far from being ~ t, on the other 
hand the data are not monotonic with d. 

One might expect that asymptotic behaviour is ob- 
served earlier for higher values of the concentration c. We 
performed therefore also 2-d simulations with c = 0.7 and 
c = 1.0. As expected, Pa and (R A ) both decreased with 
c, but the strange time dependence of (R A ) persisted. 

Finally, we show in Fig. 5 the survival probability in 
d = 3, again for c = 0.5. More precisely, we show 
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— hi Pa (t)/t. It depends quite strongly on t, but its cur- 
vature is consistent with convergence to Eq.(|l]) without 
further surprises. 

In conclusion, we have performed simulations of a trap- 
ping model in which both traps and trapped particles are 
mobile with equal mobilities. The simulations, under- 
taken in the hope of understanding the precise cross-over 
to the exactly known asymptotic behaviour, are much 
more extensive than previous simulations of this model, 
in the sense of reaching much longer times and much 
lower survival probabilities. This was possible due to a 



novel algorithm. In spite of the vastly improved numer- 
ics we have not been able to understand all details of the 
model. Some of our results are indeed quite puzzling. On 
the other hand, our methods can be possibly applied also 
to other models where absorbers are free particles undis- 
turbed by the particles which they absorb. One class of 
such models are e.g. gated reaction-diffusion systems of 
the type A+B — > with stochastically changing reaction 
rates [15[. 
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